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Abstract 

We prove polynomial-time solvability of a large class of clustering problems where a 
weighted set of items has to be partitioned into clusters with respect to some balanc¬ 
ing constraints. The data points are weighted with respect to different features and the 
clusters adhere to given lower and upper bounds on the total weight of their points with 
respect to each of these features. Further the weight-contribution of a vector to a cluster 
can depend on the cluster it is assigned to. Our interest in these types of clustering prob¬ 
lems is motivated by an application in land consolidation where the ability to perform 
this kind of balancing is crucial. 

Our framework maximizes an objective function that is convex in the summed-up 
utility of the items in each cluster. Despite hardness of convex maximization and many 
related problems, for fixed dimension and number of clusters, we are able to show that 
our clustering model is solvable in time polynomial in the number of items if the weight¬ 
balancing restrictions are defined using vectors from a hxed, hnite domain. We conclude 
our discussion with a new, efficient model and algorithm for land consolidation. 

Keywords: constrained clustering, convex maximization, integer programming, land 
consolidation 
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1. Introduction 


Partitioning a set of items while respecting some constraints is a frequent task in 
exploratory data analysis, arising in both operations research and machine learning; see 
e.g. js, 1^1 . We consider partitioning for which the sizes of the clusters are restricted with 
respect to multiple criteria. There are many applications where it is necessary to adhere 
to given bounds on the cluster sizes. 

For example, these include the modeling of polycrystals in the material sciences [l| and 
face recognition using meshes, where the original mesh is partitioned into parts of equal 
sizes to obtain an OTtimal running time for graph-theoretical methods that are applied to 
all of these parts ji]. Our interest in these types of problems comes from an application 
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in land consolidation. See for example [sl, and in particular the outreach article 
for the impact in academia and practice. 

1.1. Land consolidation 

The farmers of many agricultural and private forest regions in Bavaria and Middle 
Germany own a large number of small lots that are scattered over the whole region. The 
reasons are strict heritage laws and a frequent change of ownership. There is signihcant 
overhead driving and an unnecessarily high cost of cultivation. In such a situation, a 
land consolidation process may be initiated by the state to improve on the cost-effective 
structure of the region. Voluntary land exchanges (by means of lend-lease agreements) 
are a popular method for such a process: The existing lots are kept without changes and 
the rights of cultivation are redistributed among the farmers of the region. 

This corresponds to a combinatorial redistribution of the lots and can be modeled as a 
clustering problem where each lot is an item and each farmer is a cluster . The main goal 
is to create large connected pieces of land for each farmer. One way to do so is to represent 
the lots by their midpoints in the Euclidean plane and to use the geographical locations 
of the farmsteads of each farmer as a set of sites. Then one performs a weight-balanced 
least-squares assignment of the lots to these sites BBini. The result is a redistribution 
where the farmers’ lots he close to their farmsteads. As a positive sideeffect, many of a 
farmer’s lots are connected and can be cultivated together. 

The lots differ in several features like their size, quality of soil, shape, and attached 
subsidies, and some of these features are even different for each farmer. For example one 
farmer may be eligible for subsidies if they cultivate a given lot, while another farmer 
may not. In such a situation, the lot is more valuable for the hrst farmer. Here a 
natural constraint is that - after the redistribution - each farmer should have lots that 
(approximately) sum up to the farmer’s original total with respect to each feature. 

Of course, partitioning a weighted set of items into clusters of prescribed sizes {weight- 
balanced partitioning) is readily seen to be NP-hard, even for just two clusters and each 
item having just a single weight that is uniform for both clusters: deciding whether there 
is such a mrtition is at least as hard as the Subset Sum problem. The methods in the 
literature p, 0, Him approach this intrinsic hardness by solving least-squares problems by 
an LP relaxation and rounding. The model in B performed particularly well in practice. 

In this paper, we will present a general clustering framework that, when applied to 
land consolidation, improves on this model by dealing with its biggest shortcoming: The 
model in B (and in fact B,0[ni , too) is not able to balance the weights of clusters with 
respect to multiple features of the lots at the same time. Instead the redistribution of 
cultivation rights is done with respect to a single ‘value’ of a lot, an aggregation of all its 
properties (and this value is the same for each farmer). 

The farmers will not accept a large deviation with respect to any of the features of 
their total lots. They will only participate in the redistribution if they do not lose a 
signihcant area of land, do not lose a lot in quality of soil, and do not lose much of their 
subsidies at the same time. But the models in the literature will return ‘optimal’ solutions 
for which not even the aggregated value necessarily is within the specihed bounds (e.g. 
3%) from the original. Further, even if a farmer’s aggregated deviation is small, they may 
have received more land than before, but of much lower quality - which they will not 
accept. These are intrinsic weaknesses coming from the relaxation and rounding that are 
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performed in BBHini. In the practical implementation, this meant that a lot of the 
work still had to be done ‘by hand’. 

However, the previous methods do not use all of the favorable properties of the input 
data for agricultural regions. The data typically falls into only a hxed number of categories, 
which just comes from the way the features are measured in practice. For example, 
especially large slots (> 5 hectare) are not traded at all. Further, one does not distinguish 
between lot sizes that differ by less than a tenth of a hectare, so that one obtains a hnite 
domain of lot sizes. The quality of soil is measured with a number between 1 and 100, 
which is a hnite domain itself. But in fact in a single agricultural region it is rare to have 
more than hve different values within this range. The same happens for the subsidies 
attached to lots and other measures. With this additional assumption, we are able to 
present an algorithm that solves the problem of land consolidation exactly and efficiently. 


1.2. Our contributions 

In view of the above application, in the present paper we consider a new, generalized 
class of clustering problems where weight-balancing restrictions are dehned using vectors 
from a hxed, hnite domain. Each item specihes a vector of weights for each cluster. These 
weight vectors represent the weight the item would contribute with respect to the diherent 
features if assigned to the respective cluster. The sizes of the clusters are bounded above 
and below with respect to all features. 

Further, each item gives a vector of utility values with respect to each cluster rep¬ 
resenting the utility gained if the item is assigned to the corresponding cluster. (These 
vectors do not have to be from a hxed domain.) Each cluster ’collects’ a total utility by 
summing up the utility of all its items’ utility vectors. We then maximize an objective 
function that is convex on the sum of utility vectors of each part. 

These are the functions used for shaped partitioning B O, [IB 20|, which encompass 
many of the objective functions commonly used in data analysis tasks, such as minimal- 
variance clustering; see e.g. 0 - They are intimately related to the studies of gravity 


polytopes B) 111 - Shaped partitioning is known to be NP-hard if either the dimension of 
the utility vectors or the numbers of clusters is part of the input 1^, because it captures 
the hardness of convex maximization due to the possible existence of exponentially many 
local optima 0]. For this, the number of clusters and the dimension of the utility vectors 
are hxed in our analysis. 

Our main result is to show that this framework is polynomial-time solvable in the 
number of items, provided the weight vectors come from a hxed and hnite domain. In 
Section 2, we introduce a formal notation for these clustering models, state our main 
results, and discuss the objective functions that can be represented in our models in some 
detail. Section 3 is dedicated to the necessary proofs. They are based on polynomial-size 
and -time reductions to the maximization of a (special) convex function over a system 
of constraints with an n-fold constraint matrix. The latter can be performed efficiently, 
which comes from a combination of recent algebraic Graver bases methods [IB 000 


and geometric edge-directions and zonotope methods [17|, l20|]. We then transfer the com¬ 


plexity results to a generalized model for land consolidation, and present a new algorithm 
for it, in Section 4. 
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2. Model and results 


In the following, we partition n items {1,..., n} into p clusters tti, ..., vTp. Each item 
j has a utility matrix G Its z-th column represents a vector of d utility values 

gained if item j is assigned to cluster tt,. 

The utility of a clustering tt = (tti, ..., tt^) is /( ^ Ci,..., ^ CT), where / : —)■ 

JGtTI iSTTp 

M is a convex function presented by a comparison oracle. Note that ^ Cj is the sum of 

je-TTi 

utility vectors of the items in cluster vTj. 

Recall that a comparison oracle for / is a routine that when queried on two values y, z 
returns whether f{y) < f{z) or not. In our framework, we do not even need an explicit 
form of the function, only comparisons of its values. This makes the complexity results 
presented in the following stronger - essentially one works with as little information as 
possible. Of course, for any natural specihc convex function that arises in data analysis, 
such as a norm or squared Euclidean distances (such as the ones that appear in Section 
4), it is easy to explicitly perform these comparisons. 

The task is to hnd a clustering of maximal utility under certain balancing constraints. 
We begin with a hrst model for these balancing constraints and extend it only later on. 
First, assume that we are given a fixed, finite set O C such that each item j has a 
weight vector G O. Informally, this vector contains a weight for each cluster tt* that it 
would contribute if assigned to vr*. Further, each cluster vTj defines a total weight 6* G Z 
of the items that have to belong to it. Note that a hnite cardinality m of hi can e.g. 
be achieved by bounding the absolute values of all components in the w^: for maximal 
absolute value u, we obtain m < {2u + 1)^. 

The corresponding set of restrictions on the cluster sizes can be written as 

'^wl = bi (*e[p]). 

iSTTi 

Let n be the set of all partitions tt of {1,..., n} into p clusters. Then a full statement of 
this optimization problem (Pi) would be 

(Pi) max Ci,..., Cl) 

= h (z G [p]) 

jGTTi 

Our first main result is polynomial-time solvability of this model for fixed d and p. 

Theorem 1. Suppose d and p are fixed and suppose there is a fixed, finite set O C 
such that all G O. Then for every convex f presented by a comparison oracle, problem 
(Pi) can he solved in polynomial time. 

Note that (P) and Theorem [1] already extend the model and complexity results in 
Section 3.3 in [l3|, where only the number of items in the clusters is balanced. So essen¬ 
tially the weights of all items are 1 for each cluster. In constrast, we allow for arbitrary 
weights of the items and these weights may also differ for the different clusters. 

In many applications, instead of having exact sizes of the partition parts, one is given 
lower and upper bounds on the sizes. To extend the above program to lower and upper 
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bounds bf G Z on the total weights of the partition parts, we extend our formulation in 
several places. Formally, these constraints take the form 

jGiTi 

For our later proofs, we now rewrite the corresponding optimization problem using only 
equalities. Let us introduce slack variables sf^Z for all z G [p]. Then the new optimiza¬ 
tion problem (P 2 ) can be stated as 

(F 2 ) max /(E E Ci) 

(E^i) + E = K (ze[p]) 

(EE) - E = b- {i e [p]) 

je-Ki 

sf >0 (z e [p]) 

With some modihcations to the construction for (Pi) in the proof of Theorem [H we can 
show a similar statement for this more general class of problems. 

Theorem 2. Suppose d and p are fixed and suppose there is a fixed, finite set VL d ZP 
such that all E Q. Then for every convex f presented by a comparison oracle, problem 
(P 2 ) can be solved in polynomial time. 

Finally, we extend our model to allow for balanced weights with respect to s > 1 
different features. Instead of item j listing a vector G Z^ of weights with respect to 
the clusters, it now has a matrix G Z^^^ listing s-dimensional vectors W- of weights 
contributed to the cluster vTj it is assigned to. All of these weight matrices come from 
a hxed, hnite set hi C As before, a hnite cardinality m of Q can be achieved by 

bounding the absolute values of all components in the IFh for maximal absolute value u, 
one now obtains m < {2u + 1 )®^. 

Further, instead of bf G Z, we now use Bf" G Z®, and likewise we have to use slack 
vectors Sf G Z^ for all i E [p]. This gives us an optimization problem (P 3 ) as 

(P 3 ) max /(E Ci,..., E E) 

ieTTl iGTTp 

(Em + St = Bt 

(eV) - E = B- 

j&TTi 

St > 0 

Even for this extension, we keep polynomial-time solvability. 

Theorem 3. Suppose d and p are fixed and suppose there is a fixed, finite set hi C Z®^^ 
such that all E hi. Then for every convex f presented by a comparison oracle, problem 
(P 3 ) can be solved in polynomial time. 

We conclude this section with some general remarks on the expressive power of the ob¬ 
jective function. In general, the convex functions / : —>■ M of the form /( E E) • • •) E ^p) 

iSTTl jGTTp 


[i E \p\) 
{i E [p]) 
{i E [p]) 
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cover a wide range of objective functions commonly used in clustering. In particular, they 
directly represent the aggregation of utility values. Each of the p clusters tTj contributes 
a summed-up utility vector ^ (7/ of its items and these vectors then are aggregated to 


j€TTi 

a value in M by the convex function /. A simplemost case for / is to just add up the 
components of all these vectors possibly (scaled by a factor). If the components of all 
Cj are non-negative, such an objective function corresponds to a linear transform of the 
/i-norm. Recall that all norms are convex functions. By choosing Ig with s > 1, one values 
single components in the utility vectors of the clusters higher in relation to many equally 
large values. For example, for l^a, only the largest absolute values among all components 
counts. 

By means of our framework, one can represent even more general ways of aggregating 
the utility by turning to the so-called clustering bodies [ly : 
one II ■ II for on the parts ^ and a monotone norm 


Here one combines two norms, 
II • II* for R^ aggregating these 


j&TTi 


values. The objective function then takes the form 


f{J2 c;,..., 5^ q) = 11(115^ c;ii,..., II5^ qiDii 

je-TTi je-Kp je-Ki je-n-p 


The level set for value at most one is a convex body, a clustering body, and serves as the 
unit ball for a semi-norm. This explains convexity of the above /. 


3. Proofs 


The proofs for Theorems[T]to[3]have a common demeanor in that we exhibit polynomial¬ 
time and -size transformation of the corresponding model to a variant of shaped parti¬ 
tioning for a set of constraints dehned by an n-fold matrix. Generally speaking, we aim 
for a problem statement of the form 

max{/(C'x) : = b,l < x < u,x e (1) 

with l,u ^ Z^, C G and a convex function / : Z'^ —)■ R presented by a comparison 

oracle. Further, the n-fold matrix is derived from two matrices Ai and A 2 by the 
classical construction in the form 


A(A 


Ai ... Ai 

A 2 

A 2 


with n copies of both Ai and A 2 arranged in the depicted layout. Note that n corresponds 
to the number of items; each of them gets a set of columns with one of the building blocks 
Ai and A 2 . Clearly, N has to be a multiple of n. 

A problem of the form (1) is polynomial-time solvable if Ai and A 2 are hxed and C 
has a constant number of rows c. The following proposition sums up this important tool 
for our proofs. 

Proposition 1 f[l3l|). Let c ^ N be fixed, and let Ai, A 2 be fixed matrices. Further, let 
AA) be the n-fold matrix derived from Ai and A 2 , let l,u E Z^, let C G Z^^^, and 
/ : Z'^ —)■ R presented by a comparison oracle. 
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Then there is an algorithm that is polynomial in n and the length of binary input for 
C, I, u and b that solves 

max{/(Cx) : = b,l < x < u,x e 


In view of Proposition [H to obtain polynomial-time-solvability for our models, we will 
show that there is a polynomial-time and -size reduction to the above form in the size 
of the input for our original models. Note that a representation of the data as above is 
in stark contrast to the ‘natural’ representation in Section 2. In addition to the explicit 
prerequisites of Proposition [H namely the constant number of rows for C and the hxed 
y4i, ^ 2 , we also have to guarantee that N stays polynomial in the size of the original 
input. 

In the following, we carefully perform the necessary reductions. We begin with the 
hrst model (Pi) and Theorem [1], and then subsequently extend the constructions. 


Proof of Theorem 1. We prove the claim by transforming (Pi) to a statement in the 
form 


max{/(C'x) : = U,l < x < u,x E 


Let us begin with the constraints. Recall that m is the size of ff. Consider the two 
matrices 


w 


1 

1 


w 


m 

1 


Ai — 


wl 


Wn 




wt, 


wj, 


and 


^2 = [11... 1] e 

Further, we dehne the vectors b = {bi,, bpY and 
to write a system of equations A^^'^x = b' = f ^ V 


1 = (1 ,..., l)'^ G M" to be able 
The variables x G Pere 


correspond to decision variables in the following way: 

We have x = (x^,..., x"')^, where x^ corresponds to the j-th column-block of the n-fold 
matrix, i.e. to the j-th ’copy’ of Ai and A 2 . Each x^ takes the form x^ = (x{,..., x^)^, 
i.e. it consists of m blocks xj G X^. Note that this construction yields N = n{mp), which 
is polynomial in the input. 

Further, let us dehne lower and upper bounds on x in the form I < x < u: Choose 
lower bounds / = 0 G and dehne the upper bound vector u = ..., fo consist 

of n blocks = {u{,..., uffY G setting = 1 G Z^ if item j has weight vector tc* 
and = 0 G Z^ otherwise. 

The system 


A^^Y = 



I < X < u 


can be derived in polynomial time and is of polynomial size, as it only uses a polynomial 
number of copies of numbers from the original input (and of zeroes and ones). Let us 
discuss why it is equivalent to our original set of constraints. First, note that the lower 
and upper bounds force an integral solution x to be a 0,1-solution. 

The assignment of item j to a cluster is determined by the decision variables x^ 
corresponding to the j-th block of Ai and A 2 . The ’v42-block’ of the system tells us that 
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precisely one entry in is eqnal to 1. By the npper bounds on x, this can only be the 
case for an index which corresponds to a correct combination of cluster TTj and weight 
contribution by item j. For this, in the ’Ai-block’ of the system, the correct weight for 
Xj is added up in the equation to obtain total cluster size bi. Thus all items are assigned 
and all clusters obtain the correct total weight. 

It remains to check whether the objective function can be written in form f{Cx), 
where C has a constant number of rows. Note that 

n m 

JGtTs j = l i=l 


Thus C takes the form 


c = [(C")^... (C'y ... (CT ■ ■ ■ (c")”] e 


with m consecutive copies for each of the (C")T These are dehned as 


(cy = 


ci 


Ci 


p 


The claim now follows from the number of dp rows being constant and observing 
that Cx yields the vector (( ^ • • •, ( X] W^e satisfy all prerequisites of 

iGTTl j£lTp 

Proposition [T] and are done. □ 

Next, we prove Theorem [2] by extending the above construction. The goal is to enter 
the slack variables sf into the model while preserving the n-fold structure of the program. 
It is possible to do so by duplicating them and moving them into the building blocks of 
the matrix as follows. 


Proof of Theorem 2. Let A, b etc. refer to the construction for Theorem [TJ Their 
corresponding new counterparts are denoted by a bar, e.g. A, b. We now want a problem 
statement for (P 2 ) of the form 

mdJc{f{Cx) : A^'^y = V J < X < u, X e 


where A^') is an n-fold matrix derived by the standard construction using matrices Ai 
and A 2 . 

For this formulation, we use a vector x = (5:^,..., where the column-blocks x^ = 
{yy, y)~, x^y now also have copies of the slack variables = ((sj)’*',..., 

and yy = {{s{y,..., {sj^yy for each of the original xL Note that x G where 
N = n{2p + mp), which again is polynomial in the input, 
is set to 

r 1 


yn) ^ 


A2 


A2 






where 




E Q Ai 
0 -E 




[0 0 As] e Rix(2p+"^p) 


with E G the unit-matrix and 0 G the 0-matrix. Further b' = (b, 1)'^ with 

1 = (1,..., 1)^ G M” and b = {bf,... ,bp ,b'^,..., b~Y listing the upper and lower bounds 
on the cluster sizesji] 

Next, let us dehne lower and upper bounds T, h G for x. As before, I = 0 E 
Further, we use the upper bound vector u = (h^,..., where tA = u^)'^ G 

P n 

]^ 2 p+mp lagging with 2p entries of value = 'Yl'T! Further h* = (0,... ,0 ,m*)’^ G 

i=ij=i 

]^ 2 p+mp lagging with 2p entries 0 for all i > 1. Note that in a feasible solution no slack 
variable ever is larger than p and that p has an encoding length that is polynomial in the 
input size of (A)- 

Finally, we construct the new matrix C = (C^,..., C"') G building-blocks 

= [0 {Cy . . . [Cy] G Rc^px^p+mp^ where 0 G R‘^px 2 p jg ^ Q-matrix and there again 
are m copies of the {Cy in the blocks CK Essentially, all parts that correspond to slack 
variables are ignored for the objective function value, i.e. f{Cx) = f{Cx). 

The new system A^'^'^x = b',1 < x < u, as well as the new objective function matrix 
C, extend the original system only by introducing an additional polynomial number of 
zeroes, ones and copies of numbers from the input. The matrices Ai,A 2 still have a 
constant number of rows and columns, and C has a constant number of rows. 

It remains to discuss why the new system is equivalent to the constraints of (^ 2 )- 
Clearly, the ’x-part’ of x still has to be a 0,1-solution. These 0,1-entries have the same 
role as in the original construction, as in the A 2 -part the slack variables are only combined 
with zeroes. 

The lower bounds on the slack variables force them to be non-negative. By the upper 
bounds, only the slack variables in the hrst block of the construction can be greater than 
zero. The other slack variables are only in the system to preserve the n-fold structure 
of the matrix. Note that it suffices to allow integral values for the slack variables (recall 
X G as both the wj and the bf are integral. The rows of the Ai-blocks then guarantee 
that the cluster sizes plus the slack variables add up to the given lower and upper bounds. 
In particular, the clusters’ sizes lie in the given, bounded range. 

Thus the new system corresponds to (P 2 ) and satishes all prererequisites of Proposition 
[H in particular Ai, A 2 , and C do so. □ 


Finally, we prove Theorem |3] by extending the construction in the proof of Theorem [2] 
to allow for s different features for each item. We continue with the notation introduced 
for the proofs of both Theorem [1] and [2l 


Proof of Theorem 3. We want a problem statement of (P 3 ) of the form 


max{/(C'x) : A^^'^x = b',l < x < u,x E 


(p*) 


^In the notation of b' and in other places, we avoid double transposes when writing vectors for better 
readability when the context is clear. For example c = (a, for two column vectors a, b would be a 
column vector c. 
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where is an n-fold matrix derived by the standard constrnction for matrices Ai and 
^2- 

We nse x = {x^,... ,x"')'^, where the colnmn-blocks , {S^)~,x^)'^ all con¬ 
tain their own copies of the slack vectors (<5'-^)’'' = ,..., and (<5'-^)“ = 

{{Sl)~,..., {S0~)'^. Herewith x G for N = n{2sp + mp), which is polynomial in 
the size of the inpnt. 

For the n-fold constrnction of A^"‘\ let now 


Ai 


E 0 h ; 
0 -E A[ 


e M2.px(25p+mp) 


h; = 


Wl 


VFd 



VF2"^ 


^ j^spxmp 


V r p 


where E G is the nnit-matrix and 0 G the 0-matrix. Farther, let 


^2 = [0 0 Hs] G 

The new right-hand side vector is b' = (b, 1)^ G -yyith 1 = (1,..., 1)^ G M"" and 

b = {Bf, ..., Bp , ,..., B~)^ listing the npper and lower bonnd vectors on the clnster 

sizes. 

Again, we nse / = 0 G as lower bonnds. The npper bonnd vector u = , iEY ^ 

P n 

contains = (h,..., h, u^)'^ G M^sp+mp begins with 2p vectors \^\{ ^ 

i=l j=l 

M® (where |VF|^ refers to the vector listing the absolntes in fF/ componentwisely), and h® = 
(0,..., 0, n®)^ G ]R 2 sp-i-mp pggins with 2sp entries 0 for all i > 1. The = (n{,..., 
are dehned by setting nj = 1 G if item j has weight matrix IF* and nj = 0 G 
otherwise. 

Finally, as in the constrnction for Theorem [2l it is necessary to ignore the parts that 
correspond to slack variables for the objective fnnction valne. It is possible to do so 
by means of the matrix C = {C^,... ,C^) G R‘^p^^ that consists of n bnilding-blocks 
C® = [0 {C'Y ... (C")^] G where 0 G R‘^P^ 2 sp jg ^ Q-matrix. Then f{Cx) = 

fiCx) = fiCx). 

It remains to explain why this formnlation represents the constraints of {P 3 ). By 
dehnition of I and u, the ’x-part’ of x still is a 0,1-solntion. These 0,1-entries play the 
same role as in the original constrnction by dehnition of the A 2 . The slack variables are 
non-negative, and by the npper bonnds only the slack variables in the hrst block of the 
constrnction can be greater than zero. The components of h are snfhcently large to not 
impose a restriction. Again it snfhces to allow integral valnes for the slack variables, as 
both the IF-^ and the B^ are integral. The rows of the Ai-blocks then gnarantee that the 
clnster sizes plus the slack variables add up to the given lower and upper bounds for each 
component of the weight vectors of the clusters, which implies that all of the clusters’ 
total weights lie within the range given by the Bf. 

Thus the new system corresponds to (P 3 ) and satishes all prererequisites of Proposition 

m □ 
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4. Polynomial-time land consolidation 

Our interest in the presented clustering framework whose favorable complexity we 
studied in Sections 2 and 3 originally came from an application in land consolidation. We 
now conclude our discussion with a new model and algorithm for it. Recall the description 
of the problem in Section 1.1. 

First, let us connect the general dehnitions used for model (P 3 ) with the application. 
We here extend j^, which was a particularly useful model in practice. The items are the 
n lots which have to be divided among the p farmers tti, ..., VTp (the clusters). We begin 
with the balancing constraints. 

4-1- Weight-balancing constraints 

The lots differ in s different features, for example size, value (which is impacted by 
the quality of soil), attached subsidies, and so on. While some of these measures are 
independent of which farmer the lot is assigned to - for example the size and quality 
of soil are just hxed numbers - for others this may not be the case. The subsidies a 
farmer gets for cultivating a lot depend on several factors that differ between farmers 
- for example depending whether the farmer represents a small local family or a large 
agricultural business working in multiple regions. 

For each lot, we set up a matrix G The column vectors W- list the contribu¬ 

tion of lot j to the total of a farmer tt* with respect to the s features, if the lot is assigned 
to the farmer. The objective functions in the following will use the sizes of the lots of 
each farmer. The information on the size of lots is represented as one of the components 
in each column W- of G We give the size of lot j the explicit name uij. 

Let Bi G M® list the summed-up total features for farmer TTj. In the redistribution 
process, the farmers do not accept a large deviation with respect to any of the features. 
For an accepted change of for example 3%, one obtains B~ = 0.97 ■ Bi and Bl' = 1.03 ■ Bi. 
The accepted deviations may also be dehned differently for the different features. 

4 . 2 . Objective function 

In the following, in a generalization of the model in j^, we discuss a family of objective 
functions for land consolidation that £t with our framework. Each lot j is represented 
by its midpoint zj G Further, the farmers tt^ specify a location Vi G of their 
farmstead. The distance of a lot and a farmstead is measured by the square of their 
Euclidean distance \\vi — ZjW^. 

We begin with the classical least-squares assignment for a given, single location Vi of 
a farmstead (as a basic building block of our model). The lots are assigned to the farmers 
such that 

'^^UjWVi- Zj\\^ (/i) 

i=l 

is minimized. The sizes Uj of the lots are used as scaling factors for the distances ||nj — 
to have a fair treatment of the assignment of one large lot or of many small lots. 

It remains to explain why this objective function hts our framework, i.e. that it can 
be represented by a function / : —)■ M that is passed the p arguments of the form 

Ci and is convex on each of these sums (2ol |: 

jgTTi 
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For each item j, one can use Cj = —ojj\\vi — so that Yh ~ Yh ~ 

j&TTi je-TTi 

It then suffices to maximize the linear function / that sums up its arguments. Formally 
f{y) = l^y with 1 = (1,..., 1)'^ e Note d = 1. 

Such a least-squares assignment favors a good assignment of lots of a farmer with 
larger total size over a good assignment of a smaller farmer. Thus one may want to 
’normalize’ the different parts of this sum. We introduce this new approach as normed 
least-squares assignment. 

For this, dehne Cj = ) g and use a function / : —)■ R that sends 

p 


—LOj\\Vi—Zj\ 

UJj 

y = Cs)) E ), one obtains a 

^ ' i=l j£ni 

maximization of 


j&T!i 




UJj\\Vi — Zj\ 


je-Ki 


j£ni 


or equivalently, and more intuitively, a minimization of 

p 


1 


E xwE 


i=l .rr"' ^ j 

j£ni 


U}j\\Vi — Zj\ 


(/2 


j&TTi 


By scaling by the inverse of each farmer’s total assigned land, each farmer’s quality of 
assignment contributes equally to the hnal objective function value. 

p 

Unfortunately, f{y) = Y~^^ convex, even when restricted to a strictly positive 


i=l 


domain, so it is necessary to resort to a (provably good) approximation for / 2 . For this, 
denote the original total size of each farmers’ lots by Kj and the accepted lower and upper 
bounds on the sizes as One can then use these apriori to estimate the total size of 
lots of each farmer after the redistribution. (Compare this to the use of an approximate 
center of gravity in Uj.) This results in a minimization of 



OJi Vi 


(/s) 


which is a linear transform of a least-squares assignment and thus hts our framework, 
too. One can easily see this by using Cj = — ^||u ~ ull^- approximation error of 
optimizing in place of /2 is provably low. 

Lemma 1. Let vr be an optimal partition for and n' be optimal for / 2 . Then n is a 
(max-^)(max—)-approa;zmafffin with respect to f 2 - 

i<p i<p 

Proof. For a simple notation, we refer to the corresponding objective function values 

s-s and /2(7J‘'). As Kf < Y one has both / 2 (vr) < (max-%)/3(7r) and 

je^i 

/3(7r) < fsin') < (max^)/2(7r'), which combines to f 2 {n) < (max^) • (max^) • f 2 {'^')- 

i<p ^ i<p i<p ^ 

Thus we obtain a (max ^)(max ^)-approximation. □ 
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For an accepted upper and lower deviation of 3%, /a would yield a 1.0609-approxiniation 
error for / 2 . Note that /a is linear and is defined with d = 1, but uses the fact that 
and can differ for ii ^ i^. In the literature there are several examples for objective 
functions that use a larger d, but that only £t our framework for kJ = for all i < p. 
Examples include finding a partition of minimal variance 13|, l20| and pushing apart the 
centers of gravity of the partition parts 0,E3- Both of them can then be interpreted 
as norm-maximization over a gravity polytope or a shaped partition polytope; recall the 
hnal remarks in Section 2. 


4-3. Algorithm and Efficiency 

Theorem E] implies that for a hxed number p of farmers, and if there is a hxed set 
of vectors of lot features, the above model is solvable in polynomial time for (/i) and (/s) 
in the size of the input, in particular n. (Recall that due to the way the features of the 
lots are measured, this is not a particularly restrictive assumption.) In combination with 
Lemma [H we also obtain a polynomial-time approximation for (/ 2 ). Let us sum up this 
information. 

Theorem 4. Suppose the number of farmers p is fixed and suppose there is a fixed set 

C of weight matrices such that all G fl. Then the above model for land 

consolidation can be solved exactly in polynomial time for objective functions fi and f^. 
For /2 a [max —) -approximation can be computed in polynomial time. 

i<p i<p 

Algorithm [T] describes the method in pseudocode for objective function (/a); we use 
the notation from the proof of Theorem |3l First, the input data is used for a formulation 
(P*) of the problem in the form 

max{f{Cx) : = b',0<x<u,xe (P*) 

For better readability, we have split up the description in two parts: First, one sets up a 
problem statement in the form (P 3 ) (step 1), then this formulation is transformed to the 
form (P*) (step 2). In our practical implementation, we of course construct the matrices 
and vectors for (P*) directly. In Section 3, we showed that this can be done in polynomial 
time and that (P*) is of polynomial size in the input. 

The input to the problem is an actual distribution of farmland, which we represent by 
y G {0,with yij = 1 if lot j belongs to farmer i and ytj = 0 otherwise. In particular, 
the original distribution is feasible with respect to the given lower and upper bounds Bf. 
Thus there always is a feasible solution to (P3), and thus to (P*), as well. The proofs of 
Theorem [1] and |3] explain how to construct an initial feasible solution x^ for (P*) from y. 
This is done in step 3; the original vector is ‘expanded’ by a polynomial number of zeroes 
in the correct components. Analogously, in the hnal step 5 of the algorithm, it then is 
easy to return an optimal vector of decision variables y* G {0, from the optimal 

solution X* for (P*). 

Steps 1 to 3 provide a suitable formulation (P*) and an initial feasible solution, such 
that the problem can be solved efficiently with tools from the literature. In step 4, one 
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As a service to the reader, let us briefly explain the main ingredients of this approach. 
For a more detailed background see the survey jl^ and the books [13, 13], as well as the 
papers HQQ- The basic notion underlying these methods is the so-called Graver 


basis of the matrix Informally, this is a finite set of vectors g with some 

special properties: in particular, if a current solution x is not optimal, and / is linear, 
there is a g Q{A^A'j and an a G M such that x + a ■ g is another feasible solution and 
satisfies f{x) < f{x + a ■ g). This gives rise to an iterative augmentation scheme, which 
is performed in step 4. 

Recall that is fixed in our setting. In turn, this implies that the matrices Ai and 
A 2 are fixed (the proof of this fact was one of our main concerns in Section 3). Thus the 
Graver basis of the matrix is also fixed. The key theorem in the theory of 

n-fold integer programming states that in this situation the size of Q{A^"‘'>) is polynomial 
in n and it can be computed in time polynomial in n (see Theorem 4.4, Chapter 4, in 

0 ). 

If / is arbitrary convex, one cannot avoid computing the entire Graver basis, which is 
polynomial in n but of large degree, and use the algorithm in 13 1. However, due to (/s) 


13MI( 

tJHEI. 


being linear, it is possible to use the drastically stronger result of [1^. We briefly explain 
how this works. To show that this scheme is efficient and polynomial, it is necessary 
to check that both the number of augmentations and the computational effort for each 
augmentation is polynomial in the input. 

For the number of augmentations, it can be shown that if one performs so-called 
Graver-best augmentation steps, i.e. one always uses a best combination of g and a to 
improve on the current feasible solution, one then obtains an upper bound 0{nL) on the 
number of necessary steps, where L is the bit size of the input, that is, a bound which is 
linear in n and in L (see the proof of Lemma 3.10, Chapter 4, in 18|). 


m 


To find the Graver-best augmentation step at each iteration, the first approach, taken 


14 


was to compute the entire Graver basis, which, as said, has size polynomial in n 
but of large degree. However, in [l3| it was shown that using a sophisticated dynamic 
program in each iteration, it is possible to find a Graver-best augmentation step at that 
iteration without explicitly constructing the Graver basis, in time O(n^). This leads to 
a very fast and practically implementable total running time 0{n^L) for augmenting the 
initial feasible solution to an optimal one. 
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Algorithm 1 Polynomial-Time Land Consolidation 

Input 

• p farmers, and for each farmer i 

— accepted lower and upper bounds Bf G with respect to s features of the 
lots 

— the geographical location Vi G of a farmstead 

• n lots, and for each lot j 

— a matrix = {W(, ..., W^) G of contributions IT/ G R® with respect to 
all s features if lot j is assigned to farmer i. Uj denotes the size of lot j. The 

come from a hnite, hxed domain fl. 

— its geographical location Zj G R^ 

• Original land distribution represented by ?/ G {0, with yij = 1 if lot j belongs 

to farmer i and pij = 0 otherwise. Ki denotes the original total size of all lots of 
farmer i, and is derived from y and the uj. 

We denote by L the input size which is the total number of bits in the binary encoding 
of the data {Bf,Vi,ujj,Zj,y). The matrices are hxed and hence do not have to be 
considered for L. 

Output 

• Partition of the lots into p clusters tti, ... ,7 rp that is optimal with respect to (/a), 

i.e. an optimal solution for 


min 



UJj\\Vi 


that adheres to the accepted lower and upper bounds 

Algorithm 

1. Use slack vectors Sf G Z^, / : R^’ —)■ R with f{y) = , and matrix 

C = (C{, ..., C/) G R^^P of utility C/ = — XjW^ of lot j to describe the 

problem in the form (P 3 ) (see Section [2]) 

2. Transform the formulation (P 3 ) to the form {P*) (see the proof of Theorem [3]) 

max{/(C'a;) : = 6 ', 0 < x < h, a; G Z^} 

3. Construct a feasible solution x* for {P*) from y (see the proofs of Theorems [U and 

ED 

4. Augment the initial feasible solution x* to an optimal one x* by iteratively using 
Graver-best augmentation steps using the algorithm of [l^, in total time 0{jAL). 

5. Construct a decision variable vector y* G {0,that corresponds to the optimal 
X* and return y* (see the proofs of Theorems [H and ED 
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